Additive manufacturing audio signal processing

• Author： Chen Lequn
• Oct 2021
• A comprehensive analysis on additive manufacturing acoustic signal.
• System requirement: Python3

## Required python libraries
import pandas as pd
import numpy as np
import scipy as sp
import os
import math

# to make this notebook's output stable across runs
np.random.seed(42)

pd.plotting.register_matplotlib_converters()


import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('axes', labelsize=14)
mpl.rc('xtick', labelsize=12)
mpl.rc('ytick', labelsize=12)
%matplotlib inline
import seaborn as sns

from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report
from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
from sklearn.preprocessing import MinMaxScaler

## Audio signal processing libraries
import librosa
import librosa.display
import IPython.display as ipd

# Pytorch libraryes and torchaudio - for GPU accelerated feature extraction
import torch
import torchaudio
from torch import nn
from torch.utils.data import DataLoader
from torch.utils.data import Dataset
from torchsummary import summary

# Where to save the figures, and dataset locations
PROJECT_ROOT_DIR = "."
Audio_PATH = os.path.join(PROJECT_ROOT_DIR, "AM_audio_dataset")
ABB_LAAM_printing_path = os.path.join(Audio_PATH, "abb_laam")
kuka_wire_LAAM_printing_path = os.path.join(Audio_PATH, "kuka_wire_laam")
gas_powder_path = os.path.join(Audio_PATH, "gas_powder")
cnc_prinitng = os.path.join(Audio_PATH, "cnc_printing")
ambient_path = os.path.join(Audio_PATH, "ambient")

IMAGE_PATH = os.path.join(PROJECT_ROOT_DIR, "signal_processing_images")
os.makedirs(Audio_PATH, exist_ok=True)
os.makedirs(IMAGE_PATH, exist_ok=True)


## function for automatically save the diagram/graph into the folder 
def save_fig(fig_id, tight_layout=True, fig_extension="png", resolution=300):
    path = os.path.join(IMAGE_PATH, fig_id + "." + fig_extension)
    print("Saving figure", fig_id)
    if tight_layout:
        plt.tight_layout()
    plt.savefig(path, format=fig_extension, dpi=resolution)

# Ignore useless warnings (see SciPy issue #5998)
import warnings
warnings.filterwarnings(action="ignore", message="^internal gelsd")

Load sound from various categories and conditions

# load sounds
LAAM1, sampling_rate = librosa.load(os.path.join(ABB_LAAM_printing_path, "printing_normal_1.wav"))
LAAM2, _ = librosa.load(os.path.join(ABB_LAAM_printing_path, "printing_normal_3.wav"))
LAAM3, _ = librosa.load(os.path.join(ABB_LAAM_printing_path, "printing_normal_5.wav"))
LAAM4, _  = librosa.load(os.path.join(ABB_LAAM_printing_path, "printing4.wav"))

cnc1, _ = librosa.load(os.path.join(cnc_prinitng, "3.wav"))
cnc1, _ = librosa.load(os.path.join(cnc_prinitng, "4.wav"))
cnc1, _ = librosa.load(os.path.join(cnc_prinitng, "5.wav"))


ambient1, _ = librosa.load(os.path.join(ambient_path, "aimbient914.wav"))
ambient2, _ = librosa.load(os.path.join(ambient_path, "aimbient5.wav"))


gas1, _ = librosa.load(os.path.join(gas_powder_path, "gas2.wav"))
gas2, _ = librosa.load(os.path.join(gas_powder_path, "printing_gas_powder.wav"))

Kuka1, _ = librosa.load(os.path.join(kuka_wire_LAAM_printing_path, "9.wav"))
Kuka2, _ = librosa.load(os.path.join(kuka_wire_LAAM_printing_path, "10.wav"))
Kuka3, _ = librosa.load(os.path.join(kuka_wire_LAAM_printing_path, "8.wav"))

ipd.Audio((os.path.join(ABB_LAAM_printing_path, "printing4.wav")))

ipd.Audio((os.path.join(kuka_wire_LAAM_printing_path, "9.wav")))

ipd.Audio((os.path.join(ambient_path, "aimbient914.wav")))  

ipd.Audio((os.path.join(gas_powder_path, "gas2.wav")))  

ipd.Audio((os.path.join(cnc_prinitng, "3.wav")))

Basic information

LAAM4.shape

(661500,)

ambient1.shape

(110250,)

sampling_rate

22050

type(LAAM4)

numpy.ndarray

LAAM4_torch, sr = torchaudio.load(os.path.join(ABB_LAAM_printing_path, "printing4.wav"))

type(LAAM4_torch)

torch.Tensor

print (str(LAAM4_torch.shape[1]) + " " + str(LAAM4_torch.shape) )

1323000 torch.Size([2, 1323000])

sr

44100

SAMPLE_RATE = 22050
NUM_SAMPLES = 22050 * 10 # take 10 seconds


signal = ambient1
length_signal = signal.shape[0]
if length_signal < NUM_SAMPLES:
    num_missing_samples = NUM_SAMPLES - length_signal
    signal = np.pad(signal, [(0, num_missing_samples)], mode='constant', constant_values=0)

signal.shape

(220500,)

visualisation of raw signal

plt.figure(figsize=(15, 10))

plt.subplot(5, 1, 1)
librosa.display.waveplot(LAAM4, alpha=1)
plt.ylim((-0.05, 0.05))
plt.title("ABB LAAM printing sound")

plt.subplot(5, 1, 2)
librosa.display.waveplot(cnc1, alpha=1)
plt.ylim((-0.1, 0.1))
plt.title("CNC printing sound")

plt.subplot(5, 1, 3)
librosa.display.waveplot(ambient1, alpha=1)
plt.ylim((-0.1, 0.1))
plt.title("Ambient sound")

plt.subplot(5, 1, 4)
librosa.display.waveplot(gas1, alpha=1)
plt.ylim((-0.1, 0.1))
plt.title("gas sound")

plt.subplot(5, 1, 5)
librosa.display.waveplot(Kuka1, alpha=1)
plt.ylim((-0.1, 0.1))
plt.title("Kuka Wire based DED sound")

plt.tight_layout()


# plt.title("Time-domain visualisation of differnet sound")
plt.suptitle('Time-domain visualisation of different sound (before cutting/padding)', fontsize = 18, y=1.0005)
save_fig("Time-domain visualisation - raw data")
plt.show()

Saving figure Time-domain visualisation - raw data

Preprocessing: Cut the signal or pad the signal, unify the signal length if necessary

### preprocessing: Cut the signal, unify the signal length
SAMPLE_RATE = 22050
NUM_SAMPLES = 22050 * 10 # take 10 seconds

def cut_if_necessary(signal):
    if signal.shape[0] > NUM_SAMPLES:
        signal = signal[:NUM_SAMPLES]
    return signal

def right_pad_if_necessary(signal):
    length_signal = signal.shape[0]
    if length_signal < NUM_SAMPLES:
        num_missing_samples = NUM_SAMPLES - length_signal
        signal = np.pad(signal, [(0, num_missing_samples)], mode='constant', constant_values=0)
        # last_dim_padding = (0, num_missing_samples)
        # signal = torch.nn.functional.pad(signal, last_dim_padding)
    return signal

details regarding numpy padding method: https://stackoverflow.com/questions/35751306/python-how-to-pad-numpy-array-with-zeros/46115998

LAAM4 = cut_if_necessary(LAAM4)
cnc1 = cut_if_necessary(cnc1)
ambient1 = cut_if_necessary(ambient1)
gas1 = cut_if_necessary(gas1)
Kuka1 = cut_if_necessary(Kuka1)

LAAM4 = right_pad_if_necessary(LAAM4)
cnc1 = right_pad_if_necessary(cnc1)
ambient1 = right_pad_if_necessary(ambient1)
gas1 = right_pad_if_necessary(gas1)
Kuka1 = right_pad_if_necessary(Kuka1)

Auidio signal time-domain visualisation

plt.figure(figsize=(15, 10))

plt.subplot(5, 1, 1)
librosa.display.waveplot(LAAM4, alpha=1)
plt.ylim((-0.05, 0.05))
plt.title("ABB LAAM printing sound")

plt.subplot(5, 1, 2)
librosa.display.waveplot(cnc1, alpha=1)
plt.ylim((-0.1, 0.1))
plt.title("CNC printing sound")

plt.subplot(5, 1, 3)
librosa.display.waveplot(ambient1, alpha=1)
plt.ylim((-0.1, 0.1))
plt.title("Ambient sound")

plt.subplot(5, 1, 4)
librosa.display.waveplot(gas1, alpha=1)
plt.ylim((-0.05, 0.05))
plt.title("gas sound")

plt.subplot(5, 1, 5)
librosa.display.waveplot(Kuka1, alpha=1)
plt.ylim((-0.1, 0.1))
plt.title("Kuka Wire based DED sound")

plt.tight_layout()


# plt.title("Time-domain visualisation of differnet sound")
plt.suptitle('Time-domain visualisation of different sound (after cutting/padding)', fontsize = 18, y=1.0005)
save_fig("Time-domain visualisation of differnet sound - after cutting padding")
plt.show()

Saving figure Time-domain visualisation of differnet sound - after cutting padding

Remarks: during the experiment, the position of the microphone varies. We have not determined the optimal position for the microphone, thus, the sound amplitutes varies from sample to sample. The closer the microphone to the laser source, the sound amplitute is larger.

Amplitude envelope [time-domain feature]

• max signal amplitute within a seleted frame
• need to specify frame size and hop length(moving steps)

# sr: 22050 Hz, ROS topic rate is 10 Hz --> ROS topic segment has 2205 samples
FRAME_SIZE = 1024
HOP_LENGTH = 512

def amplitude_envelope(signal, frame_size, hop_length):
    """Calculate the amplitude envelope of a signal with a given frame size nad hop length."""
    amplitude_envelope = []
    
    # calculate amplitude envelope for each frame
    for i in range(0, len(signal), hop_length): 
        amplitude_envelope_current_frame = max(signal[i:i+frame_size]) 
        amplitude_envelope.append(amplitude_envelope_current_frame)
    
    return np.array(amplitude_envelope)   

# number of frames in amplitude envelope for 
LAAM4_ae = amplitude_envelope(LAAM4, FRAME_SIZE, HOP_LENGTH)
len(LAAM4_ae)

431

# calculate amplitude envelope cnc1, ambient1, gas1, Kuka1
cnc1_ae = amplitude_envelope(cnc1, FRAME_SIZE, HOP_LENGTH)
ambient1_ae = amplitude_envelope(ambient1, FRAME_SIZE, HOP_LENGTH)
gas1_ae = amplitude_envelope(gas1, FRAME_SIZE, HOP_LENGTH)
Kuka1_ae = amplitude_envelope(Kuka1, FRAME_SIZE, HOP_LENGTH)

len(ambient1_ae)

431

## visualisation
frames = range(len(LAAM4_ae))
t = librosa.frames_to_time(frames, hop_length=HOP_LENGTH)

# amplitude envelope is graphed in red
#----------------------------------------------
plt.figure(figsize=(15, 17))

plt.subplot(5, 1, 1)
librosa.display.waveplot(LAAM4, alpha=1)
plt.plot(t, LAAM4_ae, color="r")
plt.ylim((-0.05, 0.05))
plt.title("ABB LAAM printing sound")

plt.subplot(5, 1, 2)
librosa.display.waveplot(cnc1, alpha=1)
plt.plot(t, cnc1_ae, color="r")
plt.ylim((-0.1, 0.1))
plt.title("CNC printing sound")

plt.subplot(5, 1, 3)
librosa.display.waveplot(ambient1, alpha=1)
plt.plot(t, ambient1_ae, color="r")
plt.ylim((-0.1, 0.1))
plt.title("Ambient sound")

plt.subplot(5, 1, 4)
librosa.display.waveplot(gas1, alpha=1)
plt.plot(t, gas1_ae, color="r")
plt.ylim((-0.1, 0.1))
plt.title("gas sound")

plt.subplot(5, 1, 5)
librosa.display.waveplot(Kuka1, alpha=1)
plt.plot(t, Kuka1_ae, color="r")
plt.ylim((-0.15, 0.15))
plt.title("Kuka Wire based DED sound")

plt.tight_layout()

plt.suptitle("Time-domain acoustic amplitute envolope", fontsize = 18, y=1.0005)

save_fig("Time-domain acoustic amplitute envolope")
plt.show()

Saving figure Time-domain acoustic amplitute envolope

Root-mean-squared energy (RMSE) [time-domain feature]

rms_cnc1 = librosa.feature.rms(cnc1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
rms_ambient1 = librosa.feature.rms(ambient1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
rms_gas1 = librosa.feature.rms(gas1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
rms_Kuka1 = librosa.feature.rms(Kuka1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
rms_LAAM4 = librosa.feature.rms(LAAM4, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]

Visualisation of RMSE

frames = range(len(rms_LAAM4))
t = librosa.frames_to_time(frames, hop_length=HOP_LENGTH)

# rms energy is graphed in red
#---------------------------------------------------
plt.figure(figsize=(15, 17))

plt.subplot(5, 1, 1)
librosa.display.waveplot(LAAM4, alpha=0.7)
plt.plot(t, rms_LAAM4, color="r")
plt.ylim((-0.05, 0.05))
plt.title("ABB LAAM printing sound")

plt.subplot(5, 1, 2)
librosa.display.waveplot(cnc1, alpha=0.7)
plt.plot(t, rms_cnc1, color="r")
plt.ylim((-0.1, 0.1))
plt.title("CNC printing sound")

plt.subplot(5, 1, 3)
librosa.display.waveplot(ambient1, alpha=0.7)
plt.plot(t, rms_ambient1, color="r")
plt.ylim((-0.1, 0.1))
plt.title("Ambient sound")

plt.subplot(5, 1, 4)
librosa.display.waveplot(gas1, alpha=0.7)
plt.plot(t, rms_gas1, color="r")
plt.ylim((-0.1, 0.1))
plt.title("gas sound")

plt.subplot(5, 1, 5)
librosa.display.waveplot(Kuka1, alpha=0.7)
plt.plot(t, rms_Kuka1, color="r")
plt.ylim((-0.15, 0.15))
plt.title("Kuka Wire based DED sound")

plt.tight_layout()
plt.suptitle("Root Mean Square Energy (RMSE) (red)", fontsize = 18, y=1.0003)

save_fig("Root Mean Square Energy (RMSE) (red)")
plt.show()

Saving figure Root Mean Square Energy (RMSE) (red)

Zero-crossing rate [time-domain]

zcr_cnc1 = librosa.feature.zero_crossing_rate(cnc1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
zcr_ambient1 = librosa.feature.zero_crossing_rate(ambient1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
zcr_gas1 = librosa.feature.zero_crossing_rate(gas1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
zcr_Kuka1 = librosa.feature.zero_crossing_rate(Kuka1, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
zcr_LAAM4 = librosa.feature.zero_crossing_rate(LAAM4, frame_length=FRAME_SIZE, hop_length=HOP_LENGTH)[0]

plt.figure(figsize=(15, 5))

plt.plot(t, zcr_Kuka1, color="y", label="Kuka Wire-based AM")
plt.plot(t, zcr_cnc1, color="r", label="CNC")
plt.plot(t, zcr_ambient1, color="b", label="ambient noise")
plt.plot(t, zcr_gas1, color="g", label="gas sound")
plt.plot(t, zcr_LAAM4, color="k", label="ABB LAAM")
plt.ylim(0, 1)
plt.xlim(0, 10)
plt.legend(loc="best")
plt.title("Zero-crossing rate", fontsize = 18, pad=20)
save_fig("Zero-crossing rate")
plt.show()

Saving figure Zero-crossing rate

Discrete Fourier Transform (DFT) [Frequency-domain] - by fast fourier transform algorithm

Frequency spectrum

# derive spectrum using FT
ft = sp.fft.fft(LAAM4)
magnitude = np.absolute(ft)
frequency = np.linspace(0, sampling_rate, len(magnitude)) 

len(magnitude)

220500

sampling_rate

22050

# plot spectrum of the ABB LAAM signal
plt.figure(figsize=(18, 8))
# plt.plot(frequency[1:22050], magnitude[1:22050]) # magnitude spectrum
plt.plot(frequency[1:], magnitude[1:]) # magnitude spectrum
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.title("ABB LAAM sound Spectrum by DFT(with Nyquist frequency effect)")
save_fig("ABB LAAM sound Spectrum by DFT (showing effect of Nyquist)")
plt.show()

Saving figure ABB LAAM sound Spectrum by DFT (showing effect of Nyquist)

# plot spectrum of the ABB LAAM signal
plt.figure(figsize=(15, 5))
plt.plot(frequency[1:110250], magnitude[1:110250]) # magnitude spectrum
# plt.plot(frequency[1:], magnitude[1:]) # magnitude spectrum
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.title("ABB LAAM sound Spectrum by DFT")
save_fig("ABB LAAM sound Spectrum by DFT ")
plt.show()

Saving figure ABB LAAM sound Spectrum by DFT 

ft_LAAM4 = sp.fft.fft(LAAM4)
ft_cnc1 = sp.fft.fft(cnc1)
ft_ambient1 = sp.fft.fft(ambient1)
ft_gas1 = sp.fft.fft(gas1)
ft_Kuka1 = sp.fft.fft(Kuka1)


magnitude_LAAM4 = np.absolute(ft_LAAM4)
magnitude_cnc1 = np.absolute(ft_cnc1)
magnitude_ambient1 = np.absolute(ft_ambient1)
magnitude_Kuka1 = np.absolute(ft_Kuka1)
magnitude_gas1 = np.absolute(ft_gas1)

def plot_magnitude_spectrum(signal, sr, title, f_ratio=1):
    X = np.fft.fft(signal)
    X_mag = np.absolute(X)
    
    plt.figure(figsize=(13, 4))
    
    f = np.linspace(0, sr, len(X_mag))
    f_bins = int(len(X_mag)*f_ratio)  
    
    plt.plot(f[:f_bins], X_mag[:f_bins])
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Magnitute')
    plt.title(title)

#---------------------------------------------------


plot_magnitude_spectrum(LAAM4, sampling_rate, "ABB LAAM printing sound Spectrum", 0.1)
save_fig("ABB LAAM printing sound Spectrum")
plt.show()


plot_magnitude_spectrum(cnc1, sampling_rate, "CNC printing sound Spectrum", 0.1) # magnitude spectrum
save_fig("CNC printing sound Spectrum")
plt.show()


plot_magnitude_spectrum(ambient1, sampling_rate, "Ambient sound Spectrum", 0.1) # magnitude spectrum
save_fig("Ambient sound Spectrum")
plt.show()


plot_magnitude_spectrum(gas1, sampling_rate, "gas sound Spectrum", 0.1) # magnitude spectrum
save_fig("gas sound Spectrum")
plt.show()


plot_magnitude_spectrum(Kuka1, sampling_rate, "Kuka Wire based DED sound Spectrum", 0.1) # magnitude spectrum
save_fig("Kuka Wire based DED sound Spectrum")
plt.show()

Saving figure ABB LAAM printing sound Spectrum

Saving figure CNC printing sound Spectrum

Saving figure Ambient sound Spectrum

Saving figure gas sound Spectrum

Saving figure Kuka Wire based DED sound Spectrum

Spectrogram - Time-frequency representation

Extracting Short-Time Fourier Transform

S_LAAM4 = librosa.stft(LAAM4, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)
S_cnc1 = librosa.stft(cnc1, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)
S_ambient1 = librosa.stft(ambient1, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)
S_gas1 = librosa.stft(gas1, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)
S_Kuka1= librosa.stft(Kuka1, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)

S_LAAM4.shape

(513, 431)

type(S_LAAM4[0][0])

numpy.complex64

Calculation of spectrogram

Y_LAAM4  = np.abs(S_LAAM4) ** 2
Y_cnc1  = np.abs(S_cnc1) ** 2
Y_ambient1  = np.abs(S_ambient1) ** 2
Y_gas1  = np.abs(S_gas1) ** 2
Y_Kuka1  = np.abs(S_Kuka1) ** 2

Y_LAAM4.shape

(513, 431)

type(Y_LAAM4[0][0])

numpy.float32

Visualisation of the spectrogram (Log-Amplitute and Log frequency scale)

def plot_spectrogram(Y, sr, title, hop_length, y_axis="log"):
    plt.figure(figsize=(13, 5))
    librosa.display.specshow(Y, 
                             sr=sr, 
                             hop_length=hop_length, 
                             x_axis="time", 
                             y_axis=y_axis)
    plt.colorbar(format="%+2.f")
    plt.title(title)
    save_fig(title)

Y_log_LAAM4 = librosa.power_to_db(Y_LAAM4)
Y_log_cnc1 = librosa.power_to_db(Y_cnc1)
Y_log_ambient1 = librosa.power_to_db(Y_ambient1)
Y_log_gas1 = librosa.power_to_db(Y_gas1)
Y_log_Kuka1 = librosa.power_to_db(Y_Kuka1)

plot_spectrogram(Y_log_LAAM4, sampling_rate,"spectrogram - ABB LAAM sound", HOP_LENGTH, y_axis="log")
plot_spectrogram(Y_log_cnc1, sampling_rate, "spectrogram - CNC sound", HOP_LENGTH, y_axis="log")
plot_spectrogram(Y_log_ambient1, sampling_rate,"spectrogram - Ambient sound",  HOP_LENGTH, y_axis="log")
plot_spectrogram(Y_log_gas1, sampling_rate, "spectrogram - Gas sound", HOP_LENGTH, y_axis="log")
plot_spectrogram(Y_log_Kuka1, sampling_rate, "spectrogram - Kuka Wire based DED sound", HOP_LENGTH, y_axis="log")

Saving figure spectrogram - ABB LAAM sound
Saving figure spectrogram - CNC sound
Saving figure spectrogram - Ambient sound
Saving figure spectrogram - Gas sound
Saving figure spectrogram - Kuka Wire based DED sound

Mel-spectrogram : special time-frequency representation in mel scale

Mel filter banks - a square-shape filter function

filter_banks = librosa.filters.mel(n_fft=2048, sr=22050, n_mels=10)

filter_banks.shape

(10, 1025)

plt.figure(figsize=(10, 5))
librosa.display.specshow(filter_banks, 
                         sr=sr, 
                         x_axis="linear")
plt.colorbar(format="%+2.f")
plt.show()

mel_spectrogram_LAAM4 = librosa.feature.melspectrogram(LAAM4, sr=sampling_rate, n_fft=2048, hop_length=512, n_mels=10)
mel_spectrogram_cnc1 = librosa.feature.melspectrogram(cnc1, sr=sampling_rate, n_fft=2048, hop_length=512, n_mels=10)
mel_spectrogram_ambient1 = librosa.feature.melspectrogram(ambient1, sr=sampling_rate, n_fft=2048, hop_length=512, n_mels=10)
mel_spectrogram_gas1 = librosa.feature.melspectrogram(gas1, sr=sampling_rate, n_fft=2048, hop_length=512, n_mels=10)
mel_spectrogram_Kuka1 = librosa.feature.melspectrogram(Kuka1, sr=sampling_rate, n_fft=2048, hop_length=512, n_mels=10)

mel_spectrogram_LAAM4.shape

(10, 431)

log_mel_spectrogram_LAAM4 = librosa.power_to_db(mel_spectrogram_LAAM4) # magnitute to log scale
log_mel_spectrogram_cnc1 = librosa.power_to_db(mel_spectrogram_cnc1)
log_mel_spectrogram_ambient1 = librosa.power_to_db(mel_spectrogram_ambient1)
log_mel_spectrogram_gas1 = librosa.power_to_db(mel_spectrogram_gas1)
log_mel_spectrogram_Kuka1 = librosa.power_to_db(mel_spectrogram_Kuka1)

def plot_mel_spectrogram(log_mel_spectrogram, sr, title):
    plt.figure(figsize=(13, 3))
    librosa.display.specshow(log_mel_spectrogram, 
                             sr=sr, 
                             x_axis="time", 
                             y_axis="mel")
    plt.colorbar(format="%+2.f")
    plt.title(title)
    save_fig(title)

plot_mel_spectrogram(log_mel_spectrogram_LAAM4, sampling_rate, "Mel-spectrogram of ABB LAAM sound")
plot_mel_spectrogram(log_mel_spectrogram_cnc1, sampling_rate, "Mel-spectrogram of CNC printing sound")
plot_mel_spectrogram(log_mel_spectrogram_ambient1, sampling_rate, "Mel-spectrogram of Ambient sound")
plot_mel_spectrogram(log_mel_spectrogram_gas1, sampling_rate, "Mel-spectrogram of Gas sound")
plot_mel_spectrogram(log_mel_spectrogram_Kuka1, sampling_rate, "Mel-spectrogram of Kuka Wire-based AM sound")

Saving figure Mel-spectrogram of ABB LAAM sound
Saving figure Mel-spectrogram of CNC printing sound
Saving figure Mel-spectrogram of Ambient sound
Saving figure Mel-spectrogram of Gas sound
Saving figure Mel-spectrogram of Kuka Wire-based AM sound

Mel-Frequency Cepstrum Coefficient (MFCC)

• The most important feature used in audio signal processing
• Gathering information from mel-spectrogram time-frequency representation
• Carry the identity of the acoustic signal
• time-domain original signal -> spectrum (SDFT) -> log-spectrum -> Cepstrum (IDFT/ Discrete cosine transform)
• extracting the spectral evolope while discarding spectral detail, take first 12 - 13 coefficients
• First and second derivative (difference) of MFCCs

Extracting MFCCs

mfccs_LAAM4 = librosa.feature.mfcc(y=LAAM4, n_mfcc=13, sr=sampling_rate)
mfccs_cnc1 = librosa.feature.mfcc(y=cnc1, n_mfcc=13, sr=sampling_rate)
mfccs_ambient1 = librosa.feature.mfcc(y=ambient1, n_mfcc=13, sr=sampling_rate)
mfccs_gas1 = librosa.feature.mfcc(y=gas1, n_mfcc=13, sr=sampling_rate)
mfccs_Kuka1 = librosa.feature.mfcc(y=Kuka1, n_mfcc=13, sr=sampling_rate)

mfccs_LAAM4.shape

(13, 431)

def plot_mfcc(mfccs, sr, title):
    plt.figure(figsize=(13, 3))
    librosa.display.specshow(mfccs, 
                             x_axis="time", 
                             sr=sr)
    plt.colorbar(format="%+2.f")
    plt.title(title)
    save_fig(title)

plot_mel_spectrogram(log_mel_spectrogram_LAAM4, sampling_rate, "MFCC of ABB LAAM sound")
plot_mel_spectrogram(mfccs_cnc1, sampling_rate, "MFCC of CNC printing sound")
plot_mel_spectrogram(mfccs_ambient1, sampling_rate, "MFCC of Ambient sound")
plot_mel_spectrogram(mfccs_gas1, sampling_rate, "MFCC of Gas sound")
plot_mel_spectrogram(mfccs_Kuka1, sampling_rate, "MFCC of Kuka Wire-based AM sound")

Saving figure MFCC of ABB LAAM sound
Saving figure MFCC of CNC printing sound
Saving figure MFCC of Ambient sound
Saving figure MFCC of Gas sound
Saving figure MFCC of Kuka Wire-based AM sound

Band-Energy Ratio (frequency-domain feature)

Calculate Band Energy Ratio

def calculate_split_frequency_bin(split_frequency, sample_rate, num_frequency_bins):
    """Infer the frequency bin associated to a given split frequency."""
    
    frequency_range = sample_rate / 2
    frequency_delta_per_bin = frequency_range / num_frequency_bins
    split_frequency_bin = math.floor(split_frequency / frequency_delta_per_bin)
    return int(split_frequency_bin)

split_frequency_bin = calculate_split_frequency_bin(2000, 22050, 1025)
split_frequency_bin

185

def band_energy_ratio(spectrogram, split_frequency, sample_rate):
    """Calculate band energy ratio with a given split frequency."""
    
    split_frequency_bin = calculate_split_frequency_bin(split_frequency, sample_rate, len(spectrogram[0]))
    band_energy_ratio = []
    
    # calculate power spectrogram
    power_spectrogram = np.abs(spectrogram) ** 2
    power_spectrogram = power_spectrogram.T
    
    # calculate BER value for each frame
    for frame in power_spectrogram:
        sum_power_low_frequencies = frame[:split_frequency_bin].sum()
        sum_power_high_frequencies = frame[split_frequency_bin:].sum()
        band_energy_ratio_current_frame = sum_power_low_frequencies / sum_power_high_frequencies
        band_energy_ratio.append(band_energy_ratio_current_frame)
    
    return np.array(band_energy_ratio)

ber_LAAM4 = band_energy_ratio(S_LAAM4, 2000, sampling_rate)
ber_cnc1 = band_energy_ratio(S_cnc1, 2000, sampling_rate)
ber_ambient1 = band_energy_ratio(S_ambient1, 2000, sampling_rate)
ber_gas1 = band_energy_ratio(S_gas1, 2000, sampling_rate)
ber_Kuka1 = band_energy_ratio(S_Kuka1, 2000, sampling_rate)

d:\python37\lib\site-packages\ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in float_scalars
  from ipykernel import kernelapp as app

len(ber_LAAM4)

431

len(ber_gas1)

431

Visualise Band Energy Ratio

plt.figure(figsize=(13, 5))

plt.plot(t, ber_LAAM4, color="b",  label="ABB LAAM")
plt.plot(t, ber_cnc1, color="r", label="CNC")
plt.plot(t, ber_ambient1, color="y", label="Ambient")
plt.plot(t, ber_gas1, color="g", label="gas")
plt.plot(t, ber_Kuka1, color="k", label="Kuka Wire-based")
plt.ylim((0, 1200))
plt.legend(loc="best")
plt.title("Band Energy Ratio")
save_fig("plot of band energy ratio")
plt.show()

Saving figure plot of band energy ratio

Spectral centroid (frequency-domain feature)

• The concentration of energy

sc_LAAM4 = librosa.feature.spectral_centroid(y=LAAM4, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
sc_cnc1 = librosa.feature.spectral_centroid(y=cnc1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
sc_ambient1 = librosa.feature.spectral_centroid(y=ambient1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
sc_gas1 = librosa.feature.spectral_centroid(y=gas1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
sc_Kuka1 = librosa.feature.spectral_centroid(y=Kuka1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]

sc_Kuka1.shape

(431,)

plt.figure(figsize=(13, 5))

plt.plot(t, sc_LAAM4, color='b', label="ABB LAAM")
plt.plot(t, sc_cnc1, color='r', label="CNC printing")
plt.plot(t, sc_ambient1, color='y', label="Ambient sound")
plt.plot(t, sc_gas1, color='g', label="gas sound")
plt.plot(t, sc_Kuka1, color='k', label="Kuka Wire based AM")
plt.title("Spectral Centroid (center of gravity of sound)")
plt.xlabel("time [s]")
# plt.ylim((0, 1200))
plt.legend(loc="best")
save_fig("spectrial centriod")
plt.show()

Saving figure spectrial centriod

Spectral bandwidth - Frequency domain

• Variance to the spectral centroid

ban_LAAM4 = librosa.feature.spectral_bandwidth(y=LAAM4, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
ban_cnc1 = librosa.feature.spectral_bandwidth(y=cnc1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
ban_ambient1 = librosa.feature.spectral_bandwidth(y=ambient1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
ban_gas1 = librosa.feature.spectral_bandwidth(y=gas1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]
ban_Kuka1 = librosa.feature.spectral_bandwidth(y=Kuka1, sr=sr, n_fft=FRAME_SIZE, hop_length=HOP_LENGTH)[0]

plt.figure(figsize=(13, 5))

plt.plot(t, ban_LAAM4, color='b', label="ABB LAAM")
plt.plot(t, ban_cnc1, color='r', label="CNC printing")
plt.plot(t, ban_ambient1, color='y', label="Ambient sound")
plt.plot(t, ban_gas1, color='g', label="gas sound")
plt.plot(t, ban_Kuka1, color='k', label="Kuka Wire based AM")
plt.title("Spectral bandwidth (variance)")
plt.xlabel("time [s]")
# plt.ylim((0, 1200))
plt.legend(loc="best")
save_fig("Spectral bandwidth")
plt.show()

Saving figure Spectral bandwidth